/********************************
 * Robinson Matrix Algorithm
 *
 * Copyright (C) 2009 CNRS
 * Author : Florent AVELLANEDA, Eugen CEPOI
 * Algorithm design : Victor CHEPOI, Morgan SESTON ( http://www.lif-sud.univ-mrs.fr/%7Echepoi/robinson_chepoi_seston.pdf )
 *
 *
 * All rights reserved.
 *
 *   This file is part of Robinson Matrix Algorithm.
 *
 *   Robinson Matrix Algorithm is free software: you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation, either version 3 of the License, or
 *   (at your option) any later version.
 *
 *   Robinson Matrix Algorithm is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *   GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with Robinson Matrix Algorithm.  If not, see <http://www.gnu.org/licenses/>.
 *
 ***********************************************************************************/


#include <vector>
#include <set>
#include <cassert>
#include <algorithm>

#include "Algo.h"
#include "Graphs.h"
#include "PartialOrder.h"
#include "Util.h"

#include "AlgoGraph.h"
#include "PartialOrder.h"


using namespace std;

namespace Algo
{
    int holesPerX(std::vector<std::vector<int> > &AH, std::vector<unsigned int> &Xo, Matrix *m,
            std::vector<unsigned int> &maximalChain, double epsilon);
    bool epsilonCompatibleFor1Element(int x, Matrix *m, std::vector<unsigned int> &maximalChain,
            unsigned int leftBoundIndex, double epsilon);
    bool epsilonCompatibleFor2Elements(int x, int y, std::vector<unsigned int> &maximalChain,
            unsigned int holeX, unsigned int holeY, Matrix *m, double epsilon);
    void fillXij( std::vector<std::vector<std::set<unsigned int> > > &Xij,
            std::vector<unsigned int> &Xo, std::vector<std::vector<int> > &AH,
            std::vector<unsigned int> &maximalChain);



    std::vector<unsigned int> fittingByRobinson( Matrix &dissimilarityMatrix, Matrix &robinsonMatrix, double &outDelta, bool fast)
    {
        std::vector<double> delta = getDelta(&dissimilarityMatrix);

        dissimilarityMatrix.setFictifBounds();
        PartialOrder partialOrder(&dissimilarityMatrix);

        std::set<unsigned int> elements;
        for ( unsigned int i = 0; i < dissimilarityMatrix.getSize(); ++i )
            elements.insert(i);

        std::vector<unsigned int> totalOrder;
        std::vector<unsigned int> tmpTotalOrder;
        // Recherche dichotomique
        unsigned int i = delta.size()/2;
        unsigned int lastInd = delta.size();
        int min=0;
        int max=delta.size()-1;

        while( min <= max )
        {
            if ( !refine(elements, &dissimilarityMatrix, &partialOrder, tmpTotalOrder, delta[i], true, fast)  )
                min = i+1;
            else
            {
                outDelta = delta[i];
                max = i-1;
                lastInd = i;
                totalOrder = tmpTotalOrder;
                tmpTotalOrder.clear();

                assert( (totalOrder.size() == (dissimilarityMatrix.getSize()-2)) || !"Solution not complet" );
            }
            i = (min+max)/2;
        }

        assert( (lastInd < delta.size()) || !"Solution not found" );

        dissimilarityMatrix.removeFictifBounds();

        robinsonMatrix = TotalOrder2RobinsonMatrix(dissimilarityMatrix, totalOrder, outDelta);

        return totalOrder;
    }

    bool refine( std::set<unsigned int> &X, Matrix *m, PartialOrder *partialOrder,
            std::vector<unsigned int> &totalOrder, double epsilon, bool first, bool fast )
    {
        if( X.size()==3)
        {
            totalOrder.clear();
            totalOrder.push_back( 0 );
            return true;
        }

        if ( partialOrder->constructBeetweness() )
        {
            partialOrder->reinit(epsilon, fast);
            partialOrder->setOrderForTwoArbitraryElements( first );

            partialOrder->update(true, fast);
        }

        std::vector<unsigned int> maximalChain;
        if( !Util::calculateMaximalChain(partialOrder->getGraph(), maximalChain) )
            return false;

        if ( (maximalChain.size() == X.size()) && (maximalChain.size() > 1) )
        {
            maximalChain.pop_back();
            maximalChain.erase(maximalChain.begin());
            totalOrder = maximalChain;
            return true;
        }

        std::vector<unsigned int> Xo;
        Xo = Util::notIn( X, maximalChain);

        std::vector<std::vector<int> > AH;
        int res = 0;
        while ( res == 0  )
            res = holesPerX(AH, Xo, m, maximalChain, epsilon);

        if ( (res == -1) && (AH.size() == 0) )
            return false;

        std::vector<std::vector<std::set<unsigned int> > > Xij;
        std::vector<std::vector<std::vector<unsigned int> > > Xij_m, Xij_p;
        Xij.resize(maximalChain.size());
        Xij_m.resize(maximalChain.size());
        Xij_p.resize(maximalChain.size());

        for ( unsigned int i = 0; i < Xij.size(); ++i )
            Xij[i].resize(Xij.size());

        fillXij( Xij, Xo, AH, maximalChain);

        for ( unsigned int j = 0; j < (Xij.size()-1); ++j )
        {
            Xij_m[j].resize(Xij.size());
            Xij_p[j].resize(Xij.size());

            for ( unsigned int i = (j+1); i < (Xij[i].size()-1); ++i )
            {
                std::set<unsigned int> X_m, X_p;
                std::vector<PartialOrder*> pO_m, pO_p;
                std::vector<unsigned int> cellsOrder_m, cellsOrder_p;
                std::vector<std::set<unsigned int> > cells;

                if ( !Graphs::partitionAndSort(m, epsilon, maximalChain[j], maximalChain[j+1],
                            maximalChain[i], maximalChain[i+1], Xij[j][i], cells, X_m, X_p, pO_m, pO_p,
                            cellsOrder_m, cellsOrder_p ) )
                    return false;


                // Recursion on each cell
                for ( unsigned int k = 0; k < cellsOrder_m.size(); ++k )
                {
                    std::vector<unsigned int> totalSubOrder;
                    std::set<unsigned int> elements;
                    for ( unsigned int l = 0; l < pO_m[k]->getMatrix()->getSize(); ++l )
                        elements.insert(l);

                    if ( !refine(elements, pO_m[k]->getMatrix(), pO_m[k], totalSubOrder, epsilon, false, fast) )
                    {
                        Util::clear(pO_m);
                        Util::clear(pO_p);
                        return false;
                    }

                    std::vector<unsigned int> order = Util::transformer(totalSubOrder, pO_m[k]->getMatrix()->getColsInMatrix());
                    Xij_m[j][i].insert(Xij_m[j][i].end(), order.begin(), order.end());
                }

                for ( unsigned int k = 0; k < cellsOrder_p.size(); ++k )
                {
                    std::vector<unsigned int> totalSubOrder;
                    std::set<unsigned int> elements;
                    for ( unsigned int l = 0; l < pO_p[k]->getMatrix()->getSize(); ++l )
                        elements.insert(l);


                    if ( !refine(elements, pO_p[k]->getMatrix(), pO_p[k], totalSubOrder, epsilon, false, fast) )
                    {
                        Util::clear(pO_m);
                        Util::clear(pO_p);
                        return false;
                    }

                    std::vector<unsigned int> order = Util::transformer(totalSubOrder, pO_p[k]->getMatrix()->getColsInMatrix());
                    Xij_p[j][i].insert(Xij_p[j][i].end(), order.begin(), order.end());
                }

                Util::clear(pO_m);
                Util::clear(pO_p);
            }
        }

        // Filling the holes with Xij+ and Xij-
        for ( unsigned int i = 0; i < Xij_m.size(); ++i )
        {
            totalOrder.push_back(maximalChain[i]);
            for ( unsigned int j = 0; j < i ; ++j )
                totalOrder.insert(totalOrder.end(), Xij_p[j][i].begin(), Xij_p[j][i].end());
            for ( unsigned int j = (Xij_m.size()-1); j > i ; --j )
                totalOrder.insert(totalOrder.end(), Xij_m[i][j].begin(), Xij_m[i][j].end());
        }

        totalOrder.pop_back();
        totalOrder.erase(totalOrder.begin());

        return true;
    }

    int holesPerX(std::vector<std::vector<int> > &AH, std::vector<unsigned int> &Xo, Matrix *m, std::vector<unsigned int> &maximalChain, double epsilon)
    {
        std::vector<std::set<unsigned int> > H;
        AH.clear();
        AH.resize(m->getSize());
        H.resize( maximalChain.size()-1);

        for ( unsigned int i = 0; i < Xo.size(); ++i )
        {
            for ( unsigned int j = 0; j < (maximalChain.size()-1); ++j )
            {
                if ( epsilonCompatibleFor1Element(Xo[i], m, maximalChain, j, epsilon) )
                {
                    H[j].insert( Xo[i] );
                    AH[Xo[i]].push_back(j);
                }
            }

            if ( AH[Xo[i]].size() == 0 )
            {
                AH.clear();
                return -1;
            }
            else if ( AH[Xo[i]].size() == 1 )
            {
                maximalChain.insert(maximalChain.begin()+(*AH[Xo[i]].begin())+1, Xo[i]);
                Xo.erase(Xo.begin()+i);
                AH.clear();
                return 0;
            }
        }

        for (unsigned int i = 0; i < Xo.size(); i++)
        {
            for (unsigned int ti = 0; ti < AH[Xo[i]].size(); ti++)
            {
                for (unsigned int j = 0; j < Xo.size(); j++)
                {
                    if( AH[Xo[j]].size() == 0 )
                        continue;
                    /*
                       if(j==i)
                       continue;
                       */
                    bool compatible=false;
                    for (unsigned int tj = 0; tj < AH[Xo[j]].size(); tj++)
                    {
                        if( epsilonCompatibleFor2Elements( Xo[i], Xo[j], maximalChain, AH[Xo[i]][ti], AH[Xo[j]][tj], m, epsilon ) )
                        {
                            compatible=true;
                            break;
                        }
                    }

                    if( !compatible )
                    {
                        AH[Xo[i]].erase(AH[Xo[i]].begin()+ti);
                        ti--;

                        if ( AH[Xo[i]].size() == 0 )
                        {
                            AH.clear();
                            return -1;
                        }
                        else if ( AH[Xo[i]].size() == 1 )
                        {
                            maximalChain.insert(maximalChain.begin()+(*AH[Xo[i]].begin())+1, Xo[i]);
                            Xo.erase(Xo.begin()+i);
                            AH.clear();
                            return 0;
                        }
                        break;
                    }
                }
            }
        }

        return 1;
    }

    void fillXij( std::vector<std::vector<std::set<unsigned int> > > &Xij, std::vector<unsigned int> &Xo, std::vector<std::vector<int> > &AH, std::vector<unsigned int> &/*maximalChain*/)
    {
        for ( unsigned int i = 0; i < Xo.size(); ++i )
        {
            Xij[AH[Xo[i]][0]][AH[Xo[i]][AH[Xo[i]].size()-1]].insert(Xo[i]);
        }
    }

    bool epsilonCompatibleFor1Element(int x, Matrix *m, std::vector<unsigned int> &maximalChain, unsigned int leftBoundIndex, double epsilon)
    {
        for ( unsigned int j = 0; j <= leftBoundIndex; ++j )
            for ( unsigned int k = j+1; k <= leftBoundIndex; ++k )
                if ( (max(m->get(maximalChain[k], maximalChain[j]), m->get(x, maximalChain[k])) - 2*epsilon) > m->get(x, maximalChain[j]) )
                    return false;

        for ( unsigned int j = (leftBoundIndex+1); j < maximalChain.size() ; ++j )
            for ( unsigned int k = j+1; k < maximalChain.size(); ++k )
                if ( (max(m->get(maximalChain[k], maximalChain[j]), m->get(x, maximalChain[j])) - 2*epsilon) > m->get(x, maximalChain[k]) )
                    return false;

        if( (max( m->get( maximalChain[leftBoundIndex], x ), m->get( maximalChain[leftBoundIndex+1], x ) ) - 2*epsilon) > m->get( maximalChain[leftBoundIndex], maximalChain[leftBoundIndex+1] )  )
            return false;

        return true;
    }

    bool epsilonCompatibleFor2Elements(int x, int y, std::vector<unsigned int> &maximalChain, unsigned int holeX, unsigned int holeY, Matrix *m, double epsilon)
    {
        if(holeX > holeY)
        {
            { // swap x and y
                int tmp;

                tmp=x;
                x=y;
                y=tmp;
            }

            { // swap holeX and holeY
                unsigned int tmp;

                tmp = holeX;
                holeX=holeY;
                holeY=tmp;
            }
        }

        for ( unsigned int i = 0; i <= holeX; ++i)
            if ( (max(m->get(maximalChain[i], x), m->get(x, y)) - 2.0*epsilon) > m->get(y, maximalChain[i]) )
                return false;

        for ( unsigned int i = (holeX+1); i < holeY ; ++i)
            if ( (max(m->get(maximalChain[i], x), m->get(maximalChain[i], y)) - 2.0*epsilon) > m->get(x, y) )
                return false;

        for ( unsigned int i = (holeY+1); i < maximalChain.size() ; ++i)
            if ( (max(m->get(maximalChain[i], y), m->get(x, y)) - 2.0*epsilon) > m->get(x, maximalChain[i]) )
                return false;

        return true;
    }


    /**
     * @brief generate all delta for matrix m
     *
     * @param m matix
     *
     * @return vector with all delta
     */
    std::vector<double> getDelta(const Matrix *m)
    {
        std::set<double> tmp;
        std::vector<double> vectResult;
        std::vector<double> result;

        for ( unsigned int i = 0; i < m->getSize(); i++ )
            for ( unsigned int j = i+1; j < m->getSize(); j++ )
                tmp.insert( m->get(i, j) );

        vectResult.reserve((tmp.size()*(tmp.size()+1))/2);
        std::set<double>::iterator it2;
        for (std::set<double>::iterator it1 = tmp.begin(); it1 != tmp.end(); ++it1)
        {
            it2=it1;
            for (++it2; it2 != tmp.end(); ++it2)
                vectResult.push_back( (*it2 - *it1)/2.0 );
        }
        vectResult.push_back(0);

        std::sort (vectResult.begin(), vectResult.end());

        double last;
        for (unsigned int i = 0; i < vectResult.size();)
        {
            last = vectResult[i];

            result.push_back( last );

            do
            {
                ++i;
                if(i >= vectResult.size())
                    return result;
            } while( vectResult[i] == last );
        }

        assert( !"ERROR" );
        return result;
    }


    Matrix TotalOrder2RobinsonMatrix(const Matrix &m, std::vector<unsigned int> totalOrder, double &delta )
    {
        Matrix result = m;

        if( totalOrder.size() == 0 )
            return result;

        for (unsigned int x = 0; x < totalOrder.size(); x++)
            for (unsigned int y = x+1; y < totalOrder.size(); y++)
                result.set( x, y, m.get( totalOrder[x], totalOrder[y] ) );

        double vXY;
        unsigned int x;

        for (unsigned int fx = totalOrder.size()-1; fx > 1; fx--)
        {
            for (unsigned int y = 0; y < (totalOrder.size()-fx); y++)
            {
                x = fx+y;
                vXY = result.get( x, y );
                if( result.get( x-1, y) > vXY )
                    result.set( x-1, y, vXY);
                if( result.get( x,   y+1) > vXY )
                    result.set( x,   y+1, vXY);
            }
        }

        double maxDiff=0;
        for (unsigned int x = 0; x < totalOrder.size(); x++)
            for (unsigned int y = x+1; y < totalOrder.size(); y++)
                if( maxDiff < (m.get( totalOrder[x], totalOrder[y] ) - result.get(x, y) ))
                    maxDiff = m.get( totalOrder[x], totalOrder[y] ) - result.get(x, y);

        delta = (maxDiff/2.0);

        for (unsigned int x = 0; x < totalOrder.size(); x++)
            for (unsigned int y = x+1; y < totalOrder.size(); y++)
                result.set( x, y, result.get(x, y) + delta );

        return result;
    }
}

